##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
library(DESeq2) 

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/RSEM/CombineCounts.FilterLowExpression-MergeMutiSample.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

##########################################################################################

dat_tpm_all <- dat_tpm[,-2]

##########################################################################################
## 构建分组信息
condition <- factor(sapply(strsplit(colnames(dat_tpm_all)[-1] , "_") , "[" , 2))    
individual <- factor(sapply(strsplit(colnames(dat_tpm_all)[-1] , "_") , "[" , 1))
coldata <- data.frame(row.names = colnames(dat_tpm_all)[-1] , condition , individual)

## counts变为整数
database <- round(dat_tpm_all[,-1])
rownames(database) <- dat_tpm_all$gene_id

##########################################################################################
## 构建dds矩阵：就是利用上面的counts.txt文件（即对象database）和分组信息（即对象coldata）构建。
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition+individual)

## 利用PCA主成分分析的方法评估一下数据分布，若不同的数据较为离散则不放在一起做差异
rld <- vst(dds, blind=FALSE)
out_name <- paste0( out_path , "/plotPCA.pdf")
pdf(out_name)
plotPCA(rld, ntop=500)
dev.off()

##########################################################################################
## 标准化的矩阵给wgcna
dat_forwgcna <- varianceStabilizingTransformation(dds, blind = TRUE, fitType = "parametric")
dat_forwgcna_out <- data.frame(assay(dat_forwgcna))
dat_forwgcna_out$gene_id <- rownames(dat_forwgcna_out)

image_name <- paste0( out_path , "/CombineCounts.FilterLowExpression-MergeMutiSample.varianceStabilizingTransformation.tsv" )
write.table( dat_forwgcna_out , image_name , row.names = F , sep = "\t" , quote = F )

##########################################################################################
## deseq2两两比较的流程
deseq_workflow <- function(database_use = database_use , coldata_use = coldata_use , individual_use = individual_use , condition_use = condition_use , class1 = class1 , class2 = class2 ){

    ## 个体的差异纳入考虑
    dds_use <- DESeqDataSetFromMatrix(countData=database_use, colData=coldata_use, design=~individual_use + condition_use)
    
    ##########################################################################################
    ## 利用DESeq（）函数标准化dds矩阵；
    dds1 <- DESeq(dds_use)    # 将数据标准化，必要步骤！！！
    #resultsNames(dds1)    # 查看结果的名称。
    #dds1$condition        #默认后者的处理组比前面的对照组。
    res <- results(dds1 , contrast = c("condition_use" , class1 , class2) , cooksCutoff = FALSE )  # 必要步骤！！！
    #summary(res)          #看一下结果的概要信息，p值默认小于0.1。

    ##########################################################################################
    ## 提取差异分析结果
    # table(res$padj < 0.05)        #padj 即矫正后的P值。看看有多少差异基因满足所设的P值要求。TRUE的数值为满足要求的基因个数。
    res <- res[order(res$padj),]  #按照padj 进行升序排列
    res$class1 <- class1
    res$class2 <- class2

    return(res)

}   

##########################################################################################
## 如有多个组需要比较，建议不要将其两两分开而是一起分析，通过在results时指定contrast对象，获得两两的比较结果，
## 这样可以综合考虑所以样品中的表达计算量化因子做DESeq；
## 但是，如果你通过PCA/EDA分析发现某一组或某几组的within-group variability比其他组的大，那么还是还是两两分组分开比较吧！
## 观察到Normal比较聚集，IM、IGC、DGC比较离散，数据分布不一样，分开来进行两两比较
result_diff <- c()

for( i in 1:(length(class_type)-1)){
    class1 <- class_type[i]
    col_index1 <- grep( class1 , colnames(dat_tpm_all) , value = T )

    for( j in (i+1):length(class_type)){
        class2 <- class_type[j]
        col_index2 <- grep( class2 , colnames(dat_tpm_all) , value = T )

        database_use <- database[,c(col_index2 , col_index1)]
        condition_use <- factor(sapply(strsplit(colnames(database_use) , "_") , "[" , 2))    
        individual_use <- factor(sapply(strsplit(colnames(database_use) , "_") , "[" , 1))
        coldata_use <- data.frame(row.names = colnames(database_use) , condition_use , individual_use)

        tmp_res <- deseq_workflow(database_use = database_use , coldata_use = coldata_use , individual_use = individual_use , condition_use = condition_use , class1 = class2 , class2 = class1 )
        
        result_diff <- rbind( result_diff , tmp_res )
    }   
}

result_diff$gene_id <- rownames(result_diff)

image_name <- paste0( out_path , "/DiffGene.tsv" )
write.table( result_diff , image_name , row.names = F , sep = "\t" , quote = F )
